﻿module LCA

open System
open System.IO
open System.Collections.Generic
open MathNet.Numerics
open MathNet.Numerics.FSharp
open System.Diagnostics
open MathNet.Numerics.LinearAlgebra.Double

//set point as default decimal separator
System.Threading.Thread.CurrentThread.CurrentCulture <- System.Globalization.CultureInfo.CreateSpecificCulture("En")
System.Threading.Thread.CurrentThread.CurrentUICulture <- System.Globalization.CultureInfo.CreateSpecificCulture("En")

type ImageData = 
  {
    X : DenseVector
    Y : DenseVector
  }


let text = File.ReadAllText("semeion.data")
let dimensions = 256
let mutable optimum = false
let dataList = text.Split(' ','\n') |> Seq.toList |> List.filter(fun s -> s <> "" && s <> "\r")
let (data_vector:List<ImageData>) = ResizeArray()
let seed = Random()
let (old_gamma:List<List<float>>) = ResizeArray()
Console.WriteLine("Insert the number of clusters")
let clusters = Console.ReadLine() |> Convert.ToInt32
Console.WriteLine("Insert the tolerance to compute the algorithm")
let tolerance = Console.ReadLine() |> Convert.ToDouble
Console.WriteLine("Insert the directory of the output file")
let mutable output_directory = Console.ReadLine()
if output_directory.EndsWith("\\") |> not then
  output_directory <- output_directory + "\\"

let rec gather_data (data:string list) counter x_vector =
  let rec gather_labels (l: 'a list) steps res =
    if steps = 0 then
      (l,res)
    else
      gather_labels l.Tail (steps - 1) (res @ [l.Head])

  if data = [] then
    ()
  else
    if counter = dimensions then
      let l,y = gather_labels data 10 []
      let labels = y |> List.map(fun y1 -> y1 |> Double.Parse)
      let image_component =
        {
          X = DenseVector(x_vector |> List.toArray)
          Y = DenseVector(labels |> List.toArray)
        }
      data_vector.Add(image_component)
      gather_data l 0 []
    else
      gather_data data.Tail (counter + 1) (x_vector @ [data.Head |> Convert.ToDouble])

gather_data dataList 0 []
let n = data_vector |> Seq.length
let gamma_matrix = DenseMatrix(n, clusters)
let prob_matrix = DenseMatrix(n,clusters)

let (mu_vector:List<DenseVector>) = ResizeArray()
for i = 1 to clusters do
  let random_vector = DenseVector(n)
  for j = 1 to dimensions do
    random_vector.[j - 1] <- (seed.Next(25,75) |> float) / 100.0
  mu_vector.Add(random_vector)

let pi = DenseVector(clusters, 1.0 / (float clusters))

let pxu (x:DenseVector) (u:DenseVector) =
  let mutable prod = 1.0
  for i = 1 to dimensions do
    prod <- prod * (Math.Pow(u.[i - 1],x.[i - 1]) * Math.Pow(1.0 - u.[i - 1],1.0 - x.[i - 1]))
  prod

let gamma ng kg =
  let num() =
      pi.[kg - 1] * prob_matrix.[ng - 1,kg - 1]
  let den() =
    let mutable sum = 0.0
    for j = 1 to clusters do
      sum <- sum + (pi.[j - 1] * prob_matrix.[ng - 1,j - 1])
    sum
  let new_gamma = num() / den()
  new_gamma

let compute_nk kp =
  let mutable nk = 0.0
  for i = 1 to n do
    nk <- nk + gamma_matrix.[i - 1, kp - 1]
  nk

let pi_mean kp nk =
  nk / (float n)

let x_mean kx nk =
  let mutable sum = DenseVector(dimensions)
  for i = 1 to n do
    sum <- sum + (gamma_matrix.[i - 1, kx - 1] * data_vector.[i - 1].X)
  let new_mu = sum / nk
  for i = 1 to dimensions do
    if Math.Abs(new_mu.[i - 1] - mu_vector.[kx - 1].[i - 1]) >= tolerance then
      optimum <- false
  new_mu

let mutable iterations = 0
let timer = new Stopwatch()
Console.WriteLine("Executing LCA algorithm... It might take several seconds")
timer.Start()
while not optimum do
  iterations <- iterations + 1
  optimum <- true
  //compute p(xi|uj) and store them in a matrix for better performance
  for i = 1 to n do
    for j = 1 to clusters do
      prob_matrix.[i - 1, j - 1] <- pxu data_vector.[i - 1].X mu_vector.[j - 1]
  for i = 1 to n do
    for j = 1 to clusters do
      //E-step
      gamma_matrix.[i - 1, j - 1] <- gamma i j
  for k = 1 to clusters do
    //M-step
    let nk = compute_nk k
    pi.[k - 1] <- pi_mean k nk
    mu_vector.[k - 1] <- x_mean k nk

timer.Stop()
Console.WriteLine("Done!")

Console.WriteLine("Computing Rand index...")
let mutable a = 0
let mutable b = 0
for i = 1 to n do
  let v1 = gamma_matrix.Row(i - 1)
  let x1 = data_vector.[i - 1]
  let v1_cluster_index = v1.ToArray() |> Seq.max
  for j = 1 to n do
    let v2 = gamma_matrix.Row(j - 1)
    let x2 = data_vector.[j - 1]
    let v2_cluster_index = v2.ToArray() |> Seq.max
    if i <> j then
      if x1.Y = x2.Y && v1_cluster_index = v2_cluster_index then
        a <- a + 1
      elif x1.Y <> x2.Y && v1_cluster_index <> v2_cluster_index then
        b <- b + 1

let rand_index = 2.0 * ((float a) + (float b)) / ((float n) * ((float n) - 1.0))
Console.WriteLine("Done!")

Console.WriteLine("Creating output file...")
Directory.CreateDirectory(output_directory) |> ignore
let output_text = ResizeArray()
for i = 0 to clusters - 1 do
  output_text.Add(("u[" + i.ToString() + "] = [" + mu_vector.[i].ToString() + "]"))

File.WriteAllLines(output_directory + "lca_out_" + clusters.ToString() + ".out",(output_text |> Seq.toList) @ 
  ["Clusters = " + clusters.ToString()] @
  ["Elapsed time: " + (timer.ElapsedMilliseconds / (1000 |> int64)).ToString() + " seconds"] @
  ["Iterations: " + iterations.ToString()] @
  ["Avarage time per iteration: " + ((timer.ElapsedMilliseconds / (1000 |> int64) |> float) / (float iterations)).ToString() + " seconds" ] @
  ["Tolerance: " + tolerance.ToString()] @
  ["Rand index: " + rand_index.ToString()])
Console.WriteLine("Done!")
